Podsumowanie analizy

Dane pochodzą z Protein Data Bank, zawierają informację o ligandach. W analizie pominęliśmy kolumnę title, która zawierała powtórzone dane z innych kolumn oraz kolumnę weight_col, która składała się z wartości pustych. Dodatkowo w niektórych wierszach brakowało informacji odnośnie jednego z progów odcięcia. Predykcja ilości atomów oraz elektronów powiodła się z dość dużą trafnością co jest pokazanego w wynikach.

Wykorzystane biblioteki

library(dplyr)
library(ggplot2)
library(ggforce)
library(gganimate) #devtools::install_github('thomasp85/gganimate')
library(tidyr)
library(caret)
library(DT)
library(summarytools)
library(gifski)
library(png)
library(data.table)

Dane

Powtarzalność wyników

Aby zapewnić powtarzalność wyników ustawiamy stan losowego generatora liczb.

set.seed(123)

Wczytywanie danych

Do wczytania danych używamy funkci “fread”, aby zapewnić szybsze wczytanie danych. Kolumna title jest połączeniem kolumn pdb_code, res_name, res_id oraz chain_id, zatem możemy ją usunąć podczas wczytywania, aby uniknąć powtarzania informacji.

all_data <- fread("all_summary.csv", header = TRUE, dec=".", stringsAsFactors = FALSE) %>% select(-(blob_coverage:title))

Przetwarzanie brakujących danych

Wiersze, które w kolumnie “res_name” zawierają niepożądaną przez nas wartość zostają usunięte.

cleaned_data <- all_data %>% filter(!res_name %in% c('UNK', 'UNX', 'UNL', 'DUM', 'N', 'BLOB', 'ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'MSE', 'PHE', 'PRO', 'SEC', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'DA', 'DG', 'DT', 'DC', 'DU', 'A', 'G', 'T', 'C', 'U', 'HOH', 'H20', 'WAT'))

W zbiorze danych występuje kolumna, która równa jest NA we wszystkich wierszach.

rows_without_na_in_weight_co = filter(cleaned_data, !is.na(weight_col))
dim(rows_without_na_in_weight_co)[1]
## [1] 0

Jak widać brak wierszy, których kolumna weight_col nie ma wartości pustej, dlatego możemy ją wykluczyć z dalszej analizy. Usuwamy także wiersze zawierające wartość pustą.

cleaned_data_without_empty_col <- select(cleaned_data, -weight_col)

cleaned_data_without_empty_col <- cleaned_data_without_empty_col[complete.cases(cleaned_data_without_empty_col),]

Rozmiar danych oraz ich statystyki

Zebrane dane zawierają… 408 kolumn oraz 525666 wierszy. Dane są typu character, numeric, integer. Większość kolumn jest numeryczna. Ich podstawowe statystyki prezentują się tak:

Natomiast pozostałe kolumny są następujące:

pdb_code res_name res_id chain_id skeleton_data fo_col fc_col
Length:525666 Length:525666 Length:525666 Length:525666 Length:525666 Length:525666 Length:525666
Class :character Class :character Class :character Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character Mode :character Mode :character Mode :character

Analiza

Ograniczenie danych

Naszą analizę ograniczymy do 50 najpopularniejszych wartości kolumny res_name.

get_top_n_res_names <- function(input_data, top_count) {
input_data %>% 
  select(res_name) %>% 
  group_by(res_name) %>% 
  count() %>% 
  arrange(desc(n)) %>%
  head(top_count)
}

top_50_res_name <- get_top_n_res_names(cleaned_data_without_empty_col, 50)

data_with_most_common_res_names <- cleaned_data_without_empty_col %>% filter(res_name %in% top_50_res_name$res_name)

Rozkład jej wartości prezentuje się następująco

Korelacja

W celu sprawdzenia korelacji użyjemy korelacji Rho Spearmana, ponieważ rozkład wartości przynajmniej jednej kolumny nie jest rozkładem normalnym

correlation <- as.data.frame(as.table(cor(data_with_most_common_res_names %>% select_if(is.numeric), use="complete.obs", method="spearman")))

Usuniemy teraz korelacje kolumn samych ze sobą

correlation <- correlation %>% 
  rename(first_column = Var1, second_column = Var2, freq = Freq) %>%
  filter(first_column != second_column)

Grupujemy po pierwszej kolumnie oraz dla każdej wartości obliczamy maksymalną wartośi. Następnie wyznaczmy 10 kolumn z największą korelacją oraz filtrujemy dane do wizualizacji

top_correlated <- correlation %>% 
  group_by(first_column) %>% 
  summarise(max=max(freq, na.rm = TRUE)) %>%
  arrange(desc(max)) %>%
  head(10)

correlation <- correlation %>% filter((first_column %in% top_correlated$first_column & second_column %in% top_correlated$first_column))

Rozkład wartości atomów oraz eleketronów

Niezgodność liczby atomów

Niezgodność liczby elektronów

Rozkład wartości kolumn rozpoczynających się od “part_01”

Animacja

Wykres przedstawia ilości segmentów maski kształtu oraz maski gęstości elektronowej dla każdego progu odcięcia intensywności na podstawie danych z trzema najpopularniejszymi nazwami zasobu (res_name). Skala osi pionowej jest logarytmiczna, ponieważ różnice w wartościach są dość duże.

Regresja

Kolumny do regresji zostaną wybrane na podstawie współczynnika korelacji Rho Spearmana. Następnie dane zostaną podzielone na zbiór treningowy oraz testowy w stosunku 70:30. Model zostanie oceniony na podstawie 5-krotnej, podwójnej oceny krzyżowej. Dodatkowo dokonamy predykcji na danych testowych oraz obliczymy miary R^2 oraz RMSE.

Wyznaczanie liczby elektronów

correlation <- cor(data_with_most_common_res_names %>% select_if(is.numeric), use="complete.obs", method="spearman")

correlation <- gather(as.data.frame(correlation)['local_res_atom_non_h_electron_sum',], 'column', 'correlation_ratio') %>% 
  filter(correlation_ratio > 0.8)

electron_count_predict_data <- select(data_with_most_common_res_names, correlation$column)

indexes <- createDataPartition(electron_count_predict_data$local_res_atom_non_h_electron_sum,p=0.7, list=F)

training_data <- electron_count_predict_data[indexes,]
testing_data <- electron_count_predict_data[-indexes,]

ctrl <- trainControl(
  method = "repeatedcv",
  number = 2,
  repeats = 5)

fit <- train(local_res_atom_non_h_electron_sum ~ .,
             data = training_data,
             method = "lm",
             trControl = ctrl)
fit
## Linear Regression 
## 
## 245042 samples
##      8 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 122520, 122522, 122521, 122521, 122521, 122521, ... 
## Resampling results:
## 
##   RMSE     Rsquared  MAE      
##   0.46865  0.999972  0.1390262
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
predicted_values <- predict(fit, newdata = testing_data)


rmse_electrons <- RMSE(testing_data$local_res_atom_non_h_electron_sum, predicted_values)
r2_electrons <- R2(testing_data$local_res_atom_non_h_electron_sum, predicted_values)

Wyznaczanie liczby atomów

correlation <- cor(data_with_most_common_res_names %>% select_if(is.numeric), use="complete.obs", method="spearman")

correlation <- gather(as.data.frame(correlation)['local_res_atom_non_h_count',], 'column', 'correlation_ratio') %>% 
  filter(correlation_ratio > 0.8)

atom_count_predict_data <- select(data_with_most_common_res_names, correlation$column)

indexes <- createDataPartition(atom_count_predict_data$local_res_atom_non_h_count, p=0.7, list=F)

training_data <- atom_count_predict_data[indexes,]
testing_data <- atom_count_predict_data[-indexes,]


ctrl <- trainControl(
  method = "repeatedcv",
  number = 2,
  repeats = 5)

fit <- train(local_res_atom_non_h_count ~ .,
             data = training_data,
             method = "lm",
             trControl = ctrl)
fit
## Linear Regression 
## 
## 245042 samples
##     10 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 122521, 122521, 122522, 122520, 122521, 122521, ... 
## Resampling results:
## 
##   RMSE        Rsquared   MAE       
##   0.03612161  0.9999924  0.00517815
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
predicted_values <- predict(fit, newdata = testing_data)


rmse_atoms <- RMSE(testing_data$local_res_atom_non_h_count, predicted_values)
r2_atoms <- R2(testing_data$local_res_atom_non_h_count, predicted_values)

Wyniki

Predykcja RMSE R^2
Ilości atomów 0.0345164 0.9999931
Ilości elektronów 0.5205734 0.9999658

Klasyfikacja

top_3_res_name <- get_top_n_res_names(cleaned_data, 3)

res_name_classification_data <- data_with_most_common_res_names %>% filter(res_name %in% top_3_res_name$res_name)

res_name_classification_data <- select(res_name_classification_data, (part_00_shape_segments_count:part_02_density_Z_4_0), (resolution:FoFc_max), res_name) %>% mutate(res_name=factor(res_name))


indexes <- createDataPartition(res_name_classification_data$res_name,
                           p=0.7, list=F)


training_data <- res_name_classification_data[indexes,]
testing_data <- res_name_classification_data[-indexes,]




ctrl <- trainControl(
  method = "repeatedcv",
  number = 2,
  repeats = 5)

fit <- train(res_name ~ .,
             data = training_data,
             method = "rf",
             trControl = ctrl,
             ntree=7)
fit
## Random Forest 
## 
## 84242 samples
##   324 predictor
##     3 classes: 'EDO', 'GOL', 'SO4' 
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 42122, 42120, 42121, 42121, 42120, 42122, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##     2   0.6729873  0.4886008
##   163   0.6799530  0.5000081
##   324   0.6795162  0.4993857
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 163.
predicted_classes <- predict(fit, newdata = testing_data)


confusionMatrix(data = predicted_classes, testing_data$res_name)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   EDO   GOL   SO4
##        EDO  4445  2868  1030
##        GOL  2909  6544  1667
##        SO4  1043  1772 13824
## 
## Overall Statistics
##                                           
##                Accuracy : 0.6873          
##                  95% CI : (0.6825, 0.6921)
##     No Information Rate : 0.4576          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.5113          
##  Mcnemar's Test P-Value : 0.3107          
## 
## Statistics by Class:
## 
##                      Class: EDO Class: GOL Class: SO4
## Sensitivity              0.5294     0.5851     0.8368
## Specificity              0.8593     0.8164     0.8562
## Pos Pred Value           0.5328     0.5885     0.8308
## Neg Pred Value           0.8576     0.8143     0.8614
## Prevalence               0.2326     0.3098     0.4576
## Detection Rate           0.1231     0.1813     0.3829
## Detection Prevalence     0.2311     0.3080     0.4609
## Balanced Accuracy        0.6943     0.7007     0.8465